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Abstract —The single sensor probability hypothesis density 
(PHD) and cardinalized probability hypothesis density (CPHD) 
filters have been developed in the literature using the random 
finite set framework. The existing multisensor extensions of 
these filters have limitations such as sensor order dependence, 
numerical instahility or high computational requirements. In 
this paper we derive update equations for the multisensor 
CPHD filter. The multisensor PHD filter is derived as a special 
case. Exact implementation of the multisensor CPHD Involves 
sums over all partitions of the measurements from different 
sensors and is thus Intractable. We propose a computationally 
tractable approximation which combines a greedy measurement 
partitioning algorithm with the Gaussian mixture representation 
of the PHD. Our greedy approximation method allows the user 
to control the tradeoff between computational overhead and 
approximation accuracy. 

Index Tenns —Random finite sets, multisensor CPHD filter, 
multisensor PHD filter, multisensor multitarget tracking. 

I. Introduction 

In the multitarget tracking problem often the number of 
targets and the number of observations detected by sensors 
are unknown and time varying. Thus representing the targets 
and observations as vectors is inefficient. A finite set is a 
more suitable representation. This is the motivation for the 
random finite set framework 0 which represents target 
states and observations as realizations of random finite sets. 
The implementation of the general multitarget Bayes filter 
for random finite sets is analytically and computationally 
infeasible Q. Several approximations have been proposed 
which make suitable assumptions to derive tractable filters Q- 
0 - 

The majority of research based on random finite set theory 
has focused on single sensor multitarget tracking. The proba¬ 
bility hypothesis density (PHD) filter |[^ propagates, over time, 
the probability hypothesis density function which is defined 
over the single target state space. Improving on the PHD 
filter, the cardinalized probability hypothesis density (CPHD) 
filter 0 propagates the distribution of the number of targets 
(the cardinality) in addition to the PHD function. Various 
implementations of the PHD and CPHD filter have been 
proposed, including the Gaussian mixture implementation Q, 
Q and the sequential Monte Carlo implementation |j^. These 
algorithms have been successfully applied to the problem of 
multitarget tracking in the presence of clutter. 

A general multisensor extension of the PHD filter was 
first derived for the case of two sensors by Mahler |8), 0. 
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The filter equations were further generalized to include an 
arbitrary number of sensors by Delande et al. m- Because 
of their combinatorial nature, the exact filter update equations 
of the general multisensor PHD filter are not computationally 
tractable except for a very few simple cases. Delande et al. 
ini> derive simplifications to the filter update equations 
for the case when the fields of view of different sensors have 
limited overlap. This reduces the computational complexity 
to some extent, and a particle filter based implementation is 
presented in p^ . Jian et al. GD suggest implementing the 
general multisensor PHD filter by repeated application of the 
two sensor PHD filter Q. The implementation details for 
realizing the general multisensor PHD filter in this manner 
are not made explicit, and the reported numerical simulations 
are restricted to the case of two sensors. 

To avoid the combinatorial computational complexity of the 
general multisensor PHD filter, some approximate multisensor 
filters have been proposed in the literature. The iterated- 
corrector PHD filter 0 processes the information from dif¬ 
ferent sensors in a sequential manner. A single sensor PHD 
filter processes measurements from the first sensor. Using the 
output PHD function produced by this step as the predicted 
PHD function, another single sensor PHD filter processes 
measurements from the second sensor and so on. As a result, 
the final output strongly depends on the order in which sensors 
are processed GD . This dependence on the sensor order can be 
mitigated by employing the approximate product multisensor 
PHD and CPHD filters proposed by Mahler |15|. Although 
the final results are independent of sensor order, Ouyang and 
Ji phi have reported that Monte Carlo implementation of the 
approximate product multisensor PHD filter is unstable and 
the problem worsens as the number of sensors increases. We 
have observed a similar instability in Gaussian mixture model- 
based implementations. Ouyang and Ji fT^ have proposed a 
heuristic fix to stabilise the Monte Carlo implementation but 
it is not analytically verified. A comprehensive review of the 
different multisensor multitarget tracking algorithms based on 
random finite set theory can be found in pT] Ch. 10]. 

In this paper we derive the update equations for the general 
multisensor CPHD filter. The derivation method is similar to 
that of the general multisensor PHD filter ||^, pO) with the 
additional propagation of the cardinality distribution. The mul¬ 
tisensor CPHD filter we derive has combinatorial complexity 
and an exact implementation is computationally infeasible. 
To overcome this limitation we propose a two-step greedy 
approach based on a Gaussian mixture model implementation. 
Each step can be realized using a trellis structure constructed 
using the measurements from different sensors or measurement 
subsets for different Gaussian components. The algorithm is 
applicable to both the general multisensor CPHD and the 
general multisensor PHD filters. 
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Other trellis based algorithms have been developed for target 
tracking. For single-sensor single-target tracking, the Viterbi 
algorithm is applied over a trellis of measurements constructed 
over time in mi. Each column of the trellis is a measurement 
scan at a different time step. The Viterbi algorithm is used 
to hnd the best path in the trellis corresponding to data 
associations over time. This approach has been extended to 
multitarget tracking in m for a fixed and known number of 
targets. The nodes of the trellis correspond to different data 
association hypotheses and the transition weights are based 
on measurement likelihoods. The Viterbi algorithm was also 
applied in | |20) , in conjunction with energy based transition 
weights, to identify the iT-best non-intersecting paths over the 
measurement trellis when K targets are present. 

The form of the update equations in the general multisensor 
PHD/CPHD biters are similar to the update equations of the 
single sensor PHD/CPHD biters for extended targets ED, 
p2) . The similarity is in the sense that for extended targets 
the update equation requires partitioning of the single sensor 
measurement set which can be computationally demanding. 
Granstrom et al. pp propose a Gaussian mixture model-based 
implementation of the PHD biter for extended targets with re¬ 
duced partitioning complexity. This is done by calculating the 
Mahalanobis distance between the measurements and grouping 
together measurements which are close to each other within 
a certain threshold. Orguner et al. | |22) use a similar method 
to reduce computations in the Gaussian mixture model-based 
implementation of the CPHD biter for extended targets. 

The rest of the paper is organized as follows: Section [n| 
provides a brief overview of random bnite sets. Section [nI|for- 
mally poses the problem of multisensor multitarget tracking. In 
Section IV we summarize the prediction and update equations 
of the general multisensor CPHD biter. The derivation of 
the biter update equations is provided in the appendices. 
We present computationally tractable implementations of the 
general multisensor PHD and CPHD biters in Section |V] A 
performance comparison of the proposed biter with existing 
multisensor biters is conducted using numerical simulations in 
Section |Vl] We provide conclusions in Section |VII| 

Portions of this work are presented in a conference pa¬ 
per p4) ; the present manuscript contains detailed derivations 
and proofs which were omitted from the conference paper, 
and the present manuscript also includes a more detailed 
description and evaluation of the proposed approximation and 
implementation of the general multisensor PHD and CPHD 
biters. 


II. Background on random finite sets 
A. Random finite sets 

Random bnite sets are set-valued random variables. The 
PHD and CPHD biters are derived using notions of random 
bnite sets. This section provides a review of this background, 
introducing debnitions and notation used in the derivations 
that follow. Detailed treatments of random bnite sets and the 
related statistics in the context of multitarget tracking can be 
found in Q, E|, lIZl- 

A random bnite set is completely specibed using its prob¬ 
ability density function if it exists. The probability density 


function of a random bnite set modeling the multitarget state 
is also referred to as the multitarget density function in this 
paper. Let V be a realization of a random bnite set S with 
elements from an underlying space y, i.e. Y Q y. For the 
random bnite set S, denote its density function by /h(V). 
Let the cardinality distribution of the random bnite set S be 
Ps{n) 

=Prob(|S| = n), n= 1,2,... , (1) 

where the notation |S| denotes the cardinality of set 5. 
The probability generating function (PGF) of the cardinality 
disttibution ps{n) is debned as 

oo 

( 2 ) 

n=0 

A statistic of the random bnite set, which is used by the 
PHD and CPHD biters, is the probability hypothesis density 
(PHD) function E). For the random bnite set S debned over 
an underlying space 3^, we denote its PHD by Ds{y), y 6 3^- 
Unlike the probability density function fs{Y) which is debned 
over the space of bnite sets in 3^, the PHD function Ds{y) is 
debned over the space 3^. Instead of propagating the complete 
density function, which can be computationally challenging, 
the PHD and CPHD biters propagate the low dimensional 
PHD function over time. 

B. IIDC random finite set 

An independent and identically distributed cluster (IIDC) 
random bnite set 0 is completely specibed by its cardinality 
distribution and its spatial density function. Let S be an IIDC 
random bnite set with cardinality disttibution p-E{n) and the 
spatial density function C{y). The probability density function 
fsiY) and the PHD D^{y) of the random bnite set 5 are 
given by the relations 

/-(n = l>^|!pH(|F|) nc(y) (3) 

Ds{y) = ay)p ( 4 ) 

IJ' = E{\E\) = f^npsin). ( 5 ) 

n=0 

Samples from an IIDC random bnite set can be generated by 
brst sampling a cardinality m from its cardinality distribution 
ps{n) and then independently sampling m points from its 
spatial density function C(y)- 

The Poisson random bnite set is an example of an IIDC 
random bnite set where the cardinality distribution is assumed 
to be Poisson. The PHD biter Q models the multitarget state 
as a realization of a Poisson random bnite set and propagates 
its PHD function over time. The Poisson random bnite set 
assumption can be undesirable because the variance of the 
Poisson disttibution is equal to its mean, which implies that 
as the number of targets increases, the error in its estimation 
becomes larger. To overcome this problem the CPHD biter Q 
models the multitarget state as a realization of an IIDC random 
bnite set and propagates its PHD function and cardinality 
disttibution over time. The additional cardinality information 
allows us to more accurately model the multitarget state. 
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III. Problem formulation 
We now specify the multisensor multitarget tracking prob¬ 
lem. Let Xfc i e X he the state of the target at time 
k. In most of the tracking literature X is chosen to be 
the Euclidean space, X = K"*, where rix is the dimension 
of the single target state. If nfc > 0 targets are present 
at time k, the multitarget state can be represented by the 
finite set Xk = {x^ i,.. . Xfc Xk £ X. We assume that 
each single target state evolves according to the Markovian 
transition function fk+i\k{^k+i.i\^k,i)- New targets can arrive 
and existing targets can disappear at each time step. Let the 
survival probability of an existing target with state x at time 
k be given by the function ps«,fc(x). 

Multiple sensors make observations about the multiple tar¬ 
gets present within the monitored region. Assume that there 
are s sensors, and conditional on the multitarget state, their 
observations are independent. Measurements gathered by 
sensor j lie in the space , i.e., z^ 6 ZL Let = 
{z^ Z 2 j.,...,z:^^. be the set of measurements 

collected by the j-th sensor at time step k. The measurement 
set can be empty. We assume that each target generates at 
most one measurement per sensor at each time instant k. Each 
measurement is either associated with a target or is generated 
by the clutter process. Define Zl'^ = to be the 

collection of measurement sets gathered by all sensors at time 
k. The probability of detection of sensor j at time k is given by 
p^^(x). The function /ij fc(z|x) denotes the probability den¬ 
sity (likelihood) that sensor j makes a measurement z given 
that it detects a target with state x. Denote the probability of 
a missed detection as /j(x) = 1 - feC^)- 

The objective of the multitarget tracking problem is to form 
an estimate Xk of the multitarget state at each time step k. 
This estimate is formed using all the measurements up until 
time k obtained from all the s sensors which is denoted 
by Zl':^ = {Z|''*, Z 2 '®,..., Z^'®}. More generally, we would 
like to estimate the posterior multitarget state distribution 

fk\k{Xk\Zl;^k)- 

IV. General multisensor CPHD filter 
In this section we develop the CPHD filter equations when 
multiple sensors are present. The derivation method is similar 
to the approach used to derive the general multisensor PHD 
filter equations by Mahler |j^ and Delande et al. pO) . Since the 
CPHD filter explicitly accounts for the cardinality distribution 
of the multitarget state, the filter update equations are more 
involved. Specifically, the measurement set partitions are more 
explicitly listed when compared to the PHD filter update 
equations 0 , m- The CPHD filter also requires additional 
propagation of the cardinality distribution. 

We make the following modeling assumptions while deriv¬ 
ing the multisensor CPHD filter equations 
Assumption 1: 

a) Target birth at time A:+1 is modelled using an IIDC random 
finite set. 

b) The predicted multitarget distribution at time A: +1 is IIDC. 

c) The sensor observation processes are independent condi¬ 
tional on the multitarget state Xk+i, and the sensor clutter 
processes are IIDC. 


Before deriving the filter equations, we introduce some 
notation. Let 6^+1 (x) be the PHD function and let Pb,k+i{n) 
be the cardinality distribution of the birth process at time k. For 
the j-th sensor let Ck+i,j{z) be the clutter spatial distribution 
and let Ck+ij{t) be the PGF of the clutter cardinality distri¬ 
bution at time k + 1. Let Dk+^ki'^) denote the predicted PHD 
function and let rk+i\k{^) denote the normalized predicted 
PHD function at time k + 1 (normalized so that it integrates 
to one). Let the PGF of the predicted cardinality distribution 
Pk+i\k{'n) be denoted by Mk+i\k{t). To keep the expressions 
and derivation compact we drop the time index and write 

C/c+l,i(z) = Cj{z), Ck+l,j{t) = Ep^(x) 

Mk+l\k{t) = M{t), E (?^(x), PsvMli^) =Psv{^) 

/ij,fc+i(z|x) E /ij(z|x), rfc+i|fc(x) E r(x), mj^k+i = rrij 

P6,fe+i(n) =pf,(n), (6) 

when the time is clear from the context. Note that abbreviated 
notation is used only for convenience and the above quantities 
are in general functions of time. For functions a(y) and 6 (y), 
the notation (a, b) is defined as (a, b) = f a(y) b(y) dy. In the 
following subsections we discuss the prediction and update 
steps of the general multisensor CPHD filter. 


A. CPHD prediction step 

Since sensor information is not required in the prediction 
step, the prediction step of the CPHD filter for the multisensor 
case is the same as that for the single sensor case. Denote 
the posterior probability hypothesis density at time k as 
Dfc|fc(x) and the posterior cardinality distribution as Pk\k{n). 
The predicted probability hypothesis density function at time 
fc + 1 is given by 0 

^fe+l|fe(x) = &fc+l(x) + JPsv{yv)fk+l\k{A-^)Dk\k{^)d'N , 

(7) 


where the integral is over the complete single target state 
space. The predicted cardinality distribution at time fc + 1 is 
given by 0^0 

Pfc+i|fc(^) = 


Y,Pbh 

j=0 




l\{^k\kTPsv) {ddk\kT^ Psv) 


1-3 


{Dk\kA) 


^Pk\k (0 


( 8 ) 


where n, j and I are non-negative integers. The normalized 
predicted PHD function is given by 


/ X /X ^k+llk(^) 

(x) E rfc+i|fc(x) = 

(9) 

/^k+Hk 


where Pk+ijk = npk+ijkin). 

( 10 ) 


n=l 


B. CPHD update step 

We use the notation [l,s] to denote the set of integers 
from 1 to s. Let W £ Z^'^^ such that for all j € [l,s], 
\W\j < 1 where \W\j = |{z e IV : z € Z^^j^}|. Thus the 
subset W can have at most one measurement from each 





4 


sensor. Let W be the set of all such W. For any measurement 
subset W we can uniquely associate with it a set of pair of 
indices defined as TV = {(jJ) ’■ ^ W}. For disjoint 

subsets Wi,W 2 ,. ■ - Wn, let V = \ so that 

Wi , W 2 , • ■ ■ Wn and V partition . Think of the set Wi 
as a collection of measurements made by different sensors, all 
of which are generated by the same target and the set V as 
the collection of clutter measurements made by all the sensors. 
Let P be a partition of constructed using elements from 
the set >V and a set V, given by 


P={WuW2,...W^p\_„V}, ( 11 ) 

\p\-i 

such that IJ WiUV = Zl:-^^, (12) 

i=l 

Wi n Wj = 0, for any Wi, Wj- eP,i + j (13) 

Wi n L = 0, for any W ^P (14) 


where |P| denotes the number of elements in the partition P. 

The partition P groups the measurements in into 

disjoint subsets where each subset is either generated by a 
target (the W subsets) or generated by the clutter process (the 

V subset). Let \P\j be the number of measurements made by 
sensor j which are generated by the targets. We have; 

IP-1 

\P\j = E (15) 

i=l 

The number of measurements made by sensor j which are 
classified as clutter in the partition P is (rrij -\P\j)- Let V be 
the collection of all possible partitions P of constructed 
as above. A recursive expression for constructing the collection 

V is given in Appendix 

Denote the u*^-order derivatives of the PGFs of the clutter 
cardinality distribution and the predicted cardinality distribu¬ 
tion as 

(16) 

We use 7 to denote the probability, under the predictive 
PHD, that a target is detected by no sensor, and we thus have; 

7= f rix)flq^^{x)dx. (17) 


For concise specification of the update equations, it is useful 
to combine the terms associated with the PGF of the clutter 
cardinality distribution for a partition P. Let us define the 
quantity 

(18) 

j=i 

For a set W 6 >V and the associated index set T\y define 
the quantities 

fr{-x) I 

J del \( 






'3-{0 t*)iTw 


n 


( 19 ) 


n pXx)/iz(ziix) n 


Pw(x) = 






y'r(x)| Yl Pdi^) n 






( 20 ) 


where {j,*) indicates any pair of indices of the form 
The quantity dw can be interpreted as the ratio of the like¬ 
lihood that the measurement subset W was generated by the 
target process to the likelihood that the measurement subset W 
was generated by the clutter process. The quantity pvv(x) can 
be interpreted as the normalized pseudolikelihood contribution 
of the measurement subset W. 

The updated probability hypothesis density function 
^fc+i|fc+i (x) can be expressed as the product of the normalized 
predicted probability hypothesis density rfc+i|fe(x) at time 
k + 1 and a pseudolikelihood function. The pseudolikelihood 
function can be expressed as a linear combination of functions 
(one function for each partition P) with associated weights 
ap. The all-clutter partition P = {V} where V = has an 
associated weight ag- Define 


def P^'P 
O'O = 


°6-p \ WsP / 


^ («pM(I^|-i)( 7) n dw] 

PiV \ WiP / 


n dw 


def 

(Xp = ■ 


VKeP 


( 21 ) 


( 22 ) 


E UqM(I«I-i)( 7) n dw 

Qe-P \ WiQ j 


Note that the expression W € P only includes W 6 >V and 
does not include the component V eV of P. For the all clutter 
partition P = {V} there are no elements in the partition of the 
type W and we use the convention HwepO = 1 whenever P = 
{V}. Similary we use the convention T,w<ipi) - 9 whenever 
P={V}. 


Theorem 1. Under the conditions of Assumption 1, the 
general multisensor CPHD filter update equation for the 
probability hypothesis density is 


^/c+i|fe+i(x) 

?'/c+l|fc(x) 


S 


«o nTd(^) 

i=i 


+ E I E Pw{^) 

PeP \W^P 


(23) 


and the update equation for the posterior cardinality distribu¬ 
tion is 


Pk+l\k+l (^) 
Pk+l\k (^) 


^ ( 

PeP \ 


nl 


- 1 „| . ,m 7'~ ‘ n dw 


-IPI+1 


\ {n-\P\ +1)1 

\P\<n+l 


W^P 


^ (k;pM(I^|-i)( 7) n 

PiV \ WiP / 


(24) 


where the quantities ao, ap, pw (x) and dw ore given in 
and respectively. 

The proof of Theorem [T] is provided in Appendix It 
requires the concepts of functional derivatives, probability 
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generating functionals and the multitarget Bayes filter which 
are revised in Appendices [B| [C] and [D] respectively. The proof 
depends on an intermediate result. Lemma [T] which is proved 
in Appendix 


C. General multisensor PHD filter as a special case 


In this section we show that the general multisensor PHD 
filter can be obtained as a special case of the general multi¬ 
sensor CPHD filter when the following assumptions are made. 

Assumption 2: 

a) Target birth at time fc + 1 is modelled using a Poisson 
random finite set. 

b) The predicted multitarget distribution at time A: + 1 is 
Poisson. 

c) The sensor observation processes are independent condi¬ 
tional on the multitarget state Xk+i, and the sensor clutter 
processes are Poisson. 


Since the multitarget state distribution is modelled as Pois¬ 
son it suffices to propagate the PHD function over time. Let 
the rate of the Poisson clutter process be Xj and let Cj{z) be 
the clutter spatial distribution for the sensor. Let p,k+i\k be 
the mean predicted cardinality at time k + 1. Using the Poisson 
assumptions for the predicted multitarget distribution and the 
sensor clutter processes we have 

(25) 

cj”'(0) = AJe-^^ (26) 

Using these in ( |2l] i we have the simplification ag = p,k+i\k- 
We can also simplify the term as 




b=i 






= (g/Afe+nUT-u-E pi ) m A”" 


IPI-1 

l^kfllk 


u = l 




(28) 


■ (29) 


Since the expression (j) appears in both the 

numerator and the denominator of the term a p in the PHD up¬ 
date expression, we can ignore the portion that is independent 
of P. Hence we have 


ppM(l^l-i)(7) oc fl - 

i=i 

From ( [T9l l and ( [30l l, we can write 

ppm(i^i-i)( 7) n dw 0 

WiP WiP 


(30) 


(31) 


where dw is defined as 


dw 


del 


M/c+l|/c 



n Pd(x)/iA(z?ix) 

{i,l)^Tw 


n ^ACi(zn 

{i,l)^Tw 


n 

3-{h*)iTw 


(32) 


The PHD update equation then reduces to 

^fc+i|fc+i(x) 

?'fc+i|fc(x) 

^ ( n Y, Pw{yd) 

dk.MkWM)^Y— - - - s — (33) 

E n dw ] 

PiP \WiP / 

The above equation is equivalent to the general multisensor 
PHD update equation given in Q, pO). 


V. Implementations oe the general multisensor 
CPHD AND PHD FILTERS 

In the previous section we derived update equations for 
the general multisensor CPHD filter which propagate the 
PHD function and cardinality distribution over time. Analytic 
propagation of these quantities is difficult in general without 
imposing further conditions. In the next subsection we de¬ 
velop a Gaussian mixture-based implementation of the filter 
update equations. Although the Gaussian mixture implementa¬ 
tion is analytically tractable, it is computationally intractable. 
In Section |V-B| and |V-C| we propose greedy algorithms to 
drastically reduce computations and develop computationally 
tractable approximate implementations for the general multi¬ 
sensor CPHD and PHD filters. 


A. Gaussian mixture implementation 

We make the following assumptions to obtain closed form 
updates for equations ( |2^ and ( |24l i 
Assumption 3: 

a) The probability of detection for each sensor is constant 
throughout the single target state space; i.e., p;^(x) = 
for all X. 

b) The predicted PHD is a mixture of weighted Gaussian 
densities. 

c) The single sensor observations are linear functions of a 
single target state coiTupted by zero-mean Gaussian noise. 

d) The predicted cardinality distribution has finite support; 
i.e., there exists a positive integer ng < oo such that 
Pk+ilkin) = 0, for all n > ng. 

From the above assumptions we can express the normalized 
predicted PHD as a Gaussian mixture model 

Jk;+l\k 

r(x)= Y (34) 

2 = 1 


= 1 ; and AA(,)(x) “i’ 

is the Gaussian density function with mean and 


(0 


are non-negative 


weights satisfying 

AA v(0 


covariance matrix S 


(0 


If Hj is the observation matrix 
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for sensor j then its likelihood function can be expressed 
as /ij(z|x) = Af{z; HjXjUj). Then under the conditions 
of Assumption 3, the posterior PHD at time k + 1 can be 
expressed as a weighted mixture of Gaussian densities and 
the posterior cardinality distribution has a finite support. 

Since the probability of detection is constant we have 
7 = n^=i < 7 ^. For each partition P the quantities and 

can be easily calculated since the predicted cardinality 
distribution has finite support. The integration in the numerator 
of ( [T9l l is analytically solvable under Assumption 3 and using 
properties of Gaussian density functions Q. Hence dw can be 
analytically evaluated. From these quantities we can calculate 
ao and ap from and ( |22l l. For each measurement set 
W we can express the product r(x)pvv(x) as a sum of 
weighted Gaussian densities using the properties of Gaussian 
density functions Q. Thus from the update equation ( [23] l the 
posterior PHD can be expressed as a mixture of Gaussian 
densities. Since the predicted cardinality distribution has finite 
support, from ( |24l i, the posterior cardinality distribution also 
has finite support. Similarly, under appropriate linear Gaussian 
assumptions, the posterior PHD in ( [3^ can be expressed as a 
mixture of Gaussian densities. 

The conditions of Assumption 3 allow us to analytically 
propagate the PHD and cardinality distribution but the propa¬ 
gation is still numerically infeasible. The combinatorial nature 
of the update step can be seen from ( |2^ , ( |24l l and ( [3^ . Specif¬ 
ically, the exact implementation of the general multisensor 
CPHD and PHD filters would require evaluation of all the 
permissible partitions (i.e. all P eV) that could be constructed 
from all possible measurement subsets. The number of such 
partitions is prohibitively large and a direct implementation 
is infeasible. We now discuss an approximation of the update 
step to overcome this limitation. 

The key idea of the approximate implementation is to 
identify elements of the collection V which make a sig¬ 
nificant contribution to the update expressions. We propose 
the following two-step greedy approximation to achieve this 
within the Gaussian mixture framework. The first approxi¬ 
mation step is to select a few measurement subsets W for 
each Gaussian component. These subsets are identified by 
evaluating a score function which quantifies the likelihood that 
the subset was generated by that Gaussian component. The 
second approximation step is to greedily construct partitions 
of these subsets which are significant for the update step. The 
following subsections explain these two steps in detail. 

B. Selecting the best measurement subsets 

A measurement subset is any subset of the measurement 
set such that it contains at most one measurement per 
sensor. The total number of measurement subsets that can be 
constructed when the sensor records rrij measurements is 

S 

n + !)■ When there are many targets present and/or the 

f=i 

clutter rate is high this number can be very large. Since the size 
of the collection V depends on the number of measurement 
subsets, to develop a tractable implementation of the update 
step it is necessary to limit the number of measurement sub¬ 
sets. Instead of enumerating all possible measurement subsets. 


Sensor (j) 

12 3 4 



5 



Fig. 1. Trellis diagram for constructing measurement subsets for each 
Gaussian component. Solid blue lines represent measurement subsets retained 
after processing sensors 1 to 3. Dashed red lines correspond to extensions of 
retained measurement subsets when processing measurements from sensor 4. 


they are greedily and sequentially constructed and only a few 
are retained based on the scores associated with them. 

Consider the measurement subset W and the associated set 
Tw as defined earlier. For the Gaussian component and 
the measurement subset W we can associate a score function 
/3(*)(14/) defined as 


/3(^)(W) = 




n Pdhji^i 


zflx) 


n 

\3-(3,*)tTw , 


dx 


n ^-3 




(35) 


The above score function is obtained by splitting the dw 
term in ( [T^ for each Gaussian component. Intuitively, this 
score can be interpreted as the ratio of the likelihood that 
the measurement subset W was generated by the single target 
represented by the Gaussian component to the likelihood 
that the measurement subset W was generated by the clutter 
process. The score can be analytically calculated 

since the integral is solvable under Assumption 3 and using 
properties of Gaussian densities 0- The score is high when 
the elements of the set W truly are the measurements caused 
by the target associated with the Gaussian component. We 
use to rank measurement subsets for each Gaussian 

component and retain only a fraction of them with the highest 
scores. 

For each Gaussian component, we select the measurement 
subsets by randomly ordering the sensors and incrementally 
incorporating information from each sensor in turn. We retain 
a maximum of Wmax subsets at each step. Figure [T] provides 
a graphical representation of the algorithm in the form of 
a trellis diagram. Each column of the trellis corresponds to 
observations from one of the sensors. The sensor number is 
indicated at the top of each column. The nodes of the trellis 
correspond to the sensor observations (z}, Z 2 ,..., zf, z^,...) 
or the no detection case (z^jZ^,...). 

The process of sequential construction of measurement 
subsets can be demonstrated using an example as follows. 
The solid lines in Figure represent partial measurement 
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Gaussian component (z) 
12 3 4 




Fig. 2. Trellis diagram for constructing partitions. Solid blue lines represent 
partitions retained after processing Gaussian components 1 to 3. Red dashed 
lines correspond to partition extension when incorporating measurement 
subsets con'esponding to the 4^^ Gaussian component. 


subsets retained after processing observations from sensors 
1 to 3. Now consider the measurement subset indicated by 
the thick solid line. It corresponds to the measurement subset 
When the sensor 4 measurements are processed, 
this measurement subset is extended for each node of sensor 
4 as represented by the dashed lines. The scores 
are calculated for these new measurement subsets using the 
expression in ( [35] l but limited to only the first 4 sensors. This 
is done for each existing measurement subset in the sensor- 
measurement space and Wmax measurement subsets with high¬ 
est scores are retained and considered at the next sensor. 
Although the process of constructing measurement subsets is 
dependent on the order in which sensors are processed, we 
observe from simulations that it has no significant effect on 
filter performance. Once the subsets have been selected, the 
ordering has no further effect in the update process. 

C. Constructing partitions 

The algorithm to construct partitions from subsets is similar 
to the above algorithm used to identify the best measurement 
subsets. Since the V component of a partition is unique 
given the W components, it is sufficient to identify the W 
components to uniquely specify a partition P. A graphical 
representation of the algorithm is provided in Figure Each 
column of this trellis corresponds to the set of measurement 
subsets {Wl,W 2 , ■ ■ ■} identified by the Gaussian com¬ 
ponent. The component number (z) is indicated at the top of 
each column. The node represents the empty measurement 
subset = 0 which is always included for each component 
and it corresponds to the event that the Gaussian component 
was not detected by any of the sensors. With each valid 
partition P we associate the score dp = |~[ dw with dfz = 1. 

WiP 

We greedily identify partitions of subsets by incrementally 
incorporating measurement subsets from the different compo¬ 
nents. For example the solid lines in Figure [^correspond to the 
partitions that have been retained after processing components 


number 1 to 3. The existing partitions are expanded using the 
measurement subsets from the 4*^ component as indicated by 
the dashed lines. Some extensions are not included as they 
do not lead to a valid partition. Since the empty measurement 
subset is always included in the trellis, a partition can 
always be found. We process the Gaussian components in 
decreasing order of their associated weights. After processing 
each component, we retain a maximum of Pjnax partitions 
corresponding to the ones with highest dp. These selected 
partitions of measurement subsets are used in the update 
equations ( |2^ , ( |24| ) and ( [3^ to compute the posterior PHDs 
and cardinality distribution. In our current implementation 
we select measurement subsets and construct measurement 
partitions using only the prior PHD information. Future re¬ 
search can focus on enhancing this construction procedure by 
including the current measurements along with the prior PHD. 

For the general multisensor PHD filter a slightly more 
accurate implementation can be used and is described as 
follows. After the first approximate step of identifying mea¬ 
surement subsets for each Gaussian component, instead of the 
approximate partition construction discussed in this section, 
we can find all possible partitions from the given collection 
of measurement subsets. This problem of finding all parti¬ 
tions can be mapped to the exact cover problem p5| . An 
efficient algorithm called Dancing Links has been suggested 
by Knuth | |2^ for solving this problem. This implementation 
can be used when the number of sensors and measurement 
subsets are small. 

VI. Numerical simulations 

In this section we compare different multisensor multitarget 
tracking algorithms developed using the random finite set 
theory. Specifically we compare the following filters: iterated- 
corrector PHD (IC-PHD II)), iterated-corrector CPHD (IC- 
CPHD), general multisensor PHD (G-PHD) and the general 
multisensor CPHD (G-CPHD) filter derived in this paper. 
Models used to simulate multitarget motion and multisensor 
observations are discussed in detail in the next subsections. 
The simulated observations are used by different algorithms 
to perform multitarget tracking. All the simulations were 
conducted using MATLAB Q 


A. Target dynamics 


The single target state is a four dimensional vector x = 
[x,y,Vj;,Vy] consisting of its position coordinates x and y 
and its velocities Vx and Vy along the x-axis and y-axis 
respectively. The target state evolves according to the dis¬ 
cretized version of the continuous time nearly constant velocity 
model |27) given by 




I 0 r 0 
0 I 0 T 
0 0 10 
0 0 0 1 


Vk+l.i 


(36) 


*The MATLAB code is available at http://networks.ece.mcgill.ca/software 







the sensor is given by 


Vk+l,i ™ *^( 0 : ^7]') 



- rp3 

0 


0 


0 

rj-iS 

0 


Tfrj - 

/j-i2 

~2 

~Y 

0 

T 

0 


0 

rp'2 

n 

0 

T 


where T is the sampling period and is the intensity of the 
process noise. We simulate 100 time steps with a sampling 
period of T = Is and process noise intensity of cr^ = 0.25 to. 
Figure 3(a) shows the target tracks used in the simulations and 


Figure 3(b) shows the variation of the number of targets over 
time. All the targets originate from one of the following four 
locations (±400 to, ±400m) and targets are restricted to the 
2000 to X 2000m square region centered at the origin. Targets 
1 & 2 are present in the time range k e [1,100]; targets 3 & 4 
for k € [21,100]; targets 5 & 6 for k e [41,100]; and targets 
7 & 8 for /c 6 [61,80]. 
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(a) Target tracks 



(b) Cardinality 


Fig. 3. (a): Evolution of target tracks. The blue cross indicates origin of target 
trajectory, (b): Number of targets as function of time. 


B. Measurement model 

Measurements are collected independently by six sensors. 
When a sensor detects a target, the corresponding measure¬ 
ment consists of the position coordinates of the target cor¬ 
rupted by additive Gaussian noise. Thus if a target located at 
(x,y) is detected by a sensor, the measurement gathered by 



X 


Wx 

z = 




. y. 


. '^y . 


where Wx and Wy are independent zero-mean Gaussian noise 
terms with standard deviation and respectively. In 
our simulations we use = 10m. The probability of 

detection of each sensor is constant throughout the monitoring 
region. Five of the sensors have a hxed probability of detection 
of 0.5. The probability of detection of the sixth sensor is 
variable and is changed from 0.2 to 1 in increments of 0.1. The 
clutter measurements made by each of the sensors is Poisson 
with uniform spatial density and mean clutter rate A = 10. 


C. Filter implementation details and error metric 

All the biters model the survival probability at all times and 
at all locations as constant with psv = 0.99. The target birth 
intensity is modelled as a Gaussian mixture with four com¬ 
ponents centered at (±400, ±400,0,0), each with covariance 
matrix diag([100,100,25,25]) and weight 0.1. The target 
birth cardinality distribution is assumed Poisson with mean 
0.4. We consider two cases of sensor ordering where the sensor 
with variable probability of detection is either processed brst 
(Case 1) or last (Case 2). 

For the different multisensor biters the PHD function is 
represented by a mixture of Gaussian densities whereas the 
cardinality distribution is represented by a vector of bnite 
length which sums to one. This Gaussian mixture model ap¬ 
proximation was brst used in Q and ||^ for multitarget track¬ 
ing using single sensor PHD and CPHD biters respectively. We 
perform pruning of Gaussian components with low weights 
and merging of Gaussian components in close vicinity Q 
for computational tractability. For the iterated-corrector biters 
pruning and merging is done after processing each sensor since 
many components have negligible weight and propagating 
them has no signibcant impact on tracking accuracy. For 
the general multisensor PHD and CPHD biters pruning and 
merging is performed at the end of the update step since 
intermediate Gaussian components are not accessible. The 
general multisensor PHD and CPHD biters are implemented 
using the two-step greedy approach described in Section |V] In 
our simulations the maximum number of measurement subsets 
per Gaussian component is set to Wmax = 6 and the maximum 
number of partitions of measurement subsets is set as Pmax = 6. 
For CPHD biters, the cardinality distribution is assumed to be 
zero for n > 20. 

For the PHD biters, we estimate the number of targets by 
rounding the sum of weights of the Gaussian components to 
the nearest integer. For the CPHD biters, we estimate the 
number of targets as the peak of the posterior cardinality 
disttibution. For all the biters, the target state estimates are 
the centres of the Gaussian components with highest weights 
in the posterior PHD. To reduce the computational overhead, 
after each time step we restrict the number of Gaussian 
components to a maximum of four times the estimated number 
of targets. When the estimated number of targets is zero we 
retain a maximum of four Gaussian components. 
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The tracking performance of the different filters are com¬ 
pared using the optimal sub-pattern assignment (OSPA) error 
metric p8) . For the OSPA metric, we set the cardinality 
penalty factor c = 100 and power p = 1. The OSPA error metric 
accounts for error in estimation of target states as well as the 
error in estimation of number of targets. Given the two sets 
of estimated multitarget state and the true multitarget state, it 
finds the best permutation of the larger set which minimizes 
its distance from the smaller set and assigns a fixed penalty 
c for each cardinality error. We use the Euclidean distance 
metric and consider only the target positions while computing 
the OSPA error. 


D. Results 

The target tracks shown in Figure |3(a)| are used in all the 
simulations. The generated observation sequence is changed 
by providing a different initialization seed to the random num¬ 
ber generator. We generate 100 different observation sequences 
and report the average OSPA error obtained by running each 
multisensor filter over these 100 observation sequences. The 
probability of detection of the sensor with variable probability 


of detection is gradually increased from 0.2 to 1. Figure 4(a) 


shows the average OSPA error as the probability of detection 
is changed for the two cases. Case 1 and Case 2. The IC-PHD 
filter performs significantly worse than all the other filters. 
For the IC-PHD filter Case 1, the accuracy improves relative 
to Case 2 as the probability of detection is increased since the 
sensor with more reliable information is processed towards the 
end. 

Figure 4(b) shows a portion of Figure 4(a)| enlarged for 
clarity. We observe that for the G-PHD, IC-CPHD and G- 
CPHD filters there is very little difference between perfor¬ 
mance for Case 1 and Case 2. Thus the IC-CPHD filter 
performance does not depend significantly on the order in 
which sensors are processed. For the G-PHD and G-CPHD 
filters the order in which sensors are processed to greedily 
construct measurement subsets has little impact on the final 
filter performance. The G-CPHD filter is able to outperform 
both the G-PHD and the IC-CPHD filters and has the lowest 
average OSPA error. A box and whisker plot comparison 
of the G-PHD, IC-CPHD and G-CPHD filters is shown in 
Figure [4(^ The median OSPA error and the 25-75 percentiles 
are shown for different values of pd for the sensor with variable 
probability of detection. 

We now examine the effect of the parameters Wmax and 
Fmax, i-S-, the maximum number of measurement subsets 
and the maximum number of partitions. W^ax is varied in 
the range {1,2,4,6,8,10} and P,nax is varied in the range 
(1,2,4,6,8,10}. For this simulation we fix the probability of 
detection of all the six sensors to be 0.5. All other parameters 
of the simulation are the same as before. We do tracking using 
the same tracks as before and over 100 different observation 
sequences for each pair of (Wmaxj ^max)- 

Figure |5(a)| plots the effect of changing W^ax and Pmax 
on the average OSPA error and the average computational 
time required per time-step. Simulations were performed using 
algorithms implemented in Matlab on computers with two 





0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Pd 

(c) Box and whisker plot 


Fig. 4. (a): Average OSPA error versus the probability of detection of the 
variable sensor. The solid and dashed lines correspond to Case 1 and Case 2, 
respectively, (b): A zoomed-in version of the figure in (a) focusing on the IC- 
CPHD, G-PHD and G-CPHD filters, (c): Box and whisker plot of the OSPA 
en'or as a function of p^. Boxes indicate 25-75 interquartile range; whiskers 
extend 1.5 times the range and ‘-H’ symbols indicate outliers lying beyond the 
whiskers. 
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Xeon 4-core 2.5GHz processors and 14GB RAM. Each curve 
is obtained by fixing P^nx and changing Wmax- Dashed curves 
correspond to G-PHD filter and solid curves correspond to G- 
CPHD filters. For a given pair of (Wmax,Pmax) values both 
the filters require almost the same computational time but the 
G-CPHD filter has a lower average OSPA error compared to 
the G-PHD filter. We observe that for each curve as Wmax 
increases the average OSPA error reaches a minimum quickly 
(around Wmax = 2) and then starts rising. This is because 
as Wmax is increased the non-ideal measurement subsets also 
get involved in the construction of partitions leading to noise 
terms in the update. The computational time required grows 
approximately linearly with increase in Wmax- As Pmax is 
varied the average OSPA error saturates at around Pmax = 4 
and increasing it beyond 4 has very little impact. Increasing 
^max does not significantly raise the computational time re¬ 
quirements of the approximate G-PHD and G-CPHD filter 
implementations. 

We perform another set of simulations to study the effect of 
the number of sensors and clutter rate of the sensors on filter 
performance. In this simulation the number of sensors s is 
varied in the range {2,4,6,8,10} and the clutter rate A of the 
Poisson clutter process is varied in the range {1,5,10,15,20} 
and is same for all the sensors. We fix the probability of 
detection of all the sensors to be 0.5. The approximate greedy 
algorithm parameters are set to Wmax = 6 and Pmax = 6. All 
other parameters of the simulation are unchanged. Tracking 
is performed using the same target tracks as before and 100 
different observation sequences are generated for each pair of 
(s,A). 


Figures 5(b) and 5(c) plot the average computational time 
and the average OSPA error as the number of sensors is 
changed for different clutter rate values. Each curve is obtained 
by fixing the clutter rate and changing the number of sensors. 
Dashed curves correspond to the G-PHD filter and solid 
curves correspond to the G-CPHD filters. From Figure [5(b)| we 
observe that for approximate greedy implementations of the 
G-PHD and G-CPHD filters the computational requirements 
grow linearly with the number of sensors. As the number of 
sensors is increased the average OSPA error reduces as seen 


from Figure 5(c) The G-CPHD filter requires relatively fewer 
sensors to achieve the same accuracy as that of the G-PHD 
filter. 


E. Extension to non-linear measurement model 

In this section we extend the Gaussian mixture based filter 
implementation discussed in Section |V] to include non-linear 
measurement models using the unscented Kalman filter |29| 
approach. The unscented extensions to non-linear models 
when a single sensor is present are discussed in Q and Q 
for the PHD and CPHD filters respectively. We implement 
the unscented versions of the general multisensor PHD and 
CPHD filters by repeatedly applying the equations provided 
in 0, @. Specifically, the equations are recursively applied 
for each z € W to evaluate the score function while 

constructing the measurement subsets. 

As an example, we consider the setup described in (30| 
based on at-sea experiments. Two targets are present within 
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(a) Solid lines: G-CPHD; Dashed lines: G-PHD 




(c) Solid lines: G-CPHD; Dashed lines: G-PHD 

Fig. 5. (a): Average OSPA error versus computational time per time-step 

obtained by changing VFmax in the range {1,2,4,6,8,10} and Pmax in the 
range (1, 2,4,6,8,10}. Blue dashed curves correspond to G-PHD filter and 
red solid curves correspond to G-CPHD filter, (b): Computation time required 
as a function of number of sensors, (c): Average OSPA eiTor as a function of 
number of sensors. 
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the monitoring region and portions of their tracks are shown 
in Figure 6(b) The target state x = [x,t/] consists of its 


coordinates in the x-y plane and the hlters model the motion 
of individual targets using a random walk model given by 
= Xfc^i + r]k+i^i where the process noise rik+i,i is zero- 
mean Gaussian with covariance matrix = cr^diag(l, 1). In 
our simulations we set = 0.24km. Although we consider 
linear target dynamics in this paper, the unscented approach 
can be easily extended to include non-linear target dynamics 
as well. 

The targets are monitored using acoustic sensors which 
collect bearings (angle) measurements. If sensor j is present 
at location [x-^ , j/-^] and a target detected by the sensor has 
coordinates [x, y] then the measurement made by this sensor 
is given by 


z = arctan 


(^) 

\ X - x-t ) 


+ w 


(38) 


where the measurement noise w is zero mean Gaussian with 
standard deviation and ‘arctan’ denotes the four-quadrant 
inverse tangent function. The sensor locations are assumed 
to be known. The measurements z are in the range [0,360) 
degrees. Along with the target related measurements the 
sensors also record clutter measurements not associated with 
any target. Five sensors (which slowly drift over time) gather 
measurements and their approximate locations are indicated in 
Figure 6(b) To demonstrate the feasibility of the proposed al¬ 
gorithms for non-linear measurement models, in this paper we 
consider sensor deployments and target trajectories based on 
at-sea experiments. The sensor measurements themselves are 
simulated to avoid the issue of measurement model mismatch. 
All sensors are assumed to have the same CTu, and we vary 
ayj in the range {1,2,3,4} (degrees) in our simulations. The 
probability of detection of each sensor is uniform throughout 
the monitoring region and is the same for all the sensors. The 
probability of detection of the sensors is changed from 0.7 to 
0.95 in increments of 0.05. The clutter measurements made by 
each of the sensors is a Poisson random hnite set with uniform 
density in [0,360) and mean clutter rate A = 5. 

The general multisensor PHD and the general multisensor 
CPHD hlters are used to perform tracking in this setup. Most 
of the implementation details are the same as discussed in 
Section |VI-C The target birth intensity is modeled as a two 
component Gaussian mixture with components centered at true 
location of the targets at time fc = I and each having covariance 
matrix diag([0.65,0.79]) and weight O.I. The target birth 
cardinality distribution is assumed to be Poisson with mean 
0.2. We set IFniax = 6 and Pmax = 6. While calculating the 
OSPA error we use the cardinality penalty factor of c = 2 and 
power p = I. To demonstrate the feasibility of the proposed 
algorithm, in our current implementation a-priori information 
about initial target locations is assumed to be known. In a 
more practical scenario where this information is unavailable 
the Gaussian components (for birth) can be initialized based 
on sensor measurements. 

The average OSPA error obtained by running the algorithms 
over 100 different observation sequences are shown in Fig¬ 
ure 6(a) Each curve is obtained by varying the probability of 


detection of the sensors from 0.7 to 0.95. As the probability 
of detection increases there is gradual decrease in the average 
OSPA error. Different curves correspond to different values of 
measurement noise standard deviation aw As a^ is increased 
the average OSPA error increases as expected. For each a^, 
the G-CPHD hlter performs better than the G-PHD hlter. 
Estimated target locations obtained by the general multisensor 
CPHD hlter are shown in Eigure 6(b) when a^j = 2 and 
Pd = 0.9. The tracks are obtained by joining the closest 
estimates across time. 




(b) Solid lines: True tracks; Dashed lines: Estimated tracks 


Fig. 6. (a): Average OSPA error Vs probability of detection of individual sen¬ 
sors. Different plots obtained by changing Cu, in the range {1,2,3,4}. Blue 
dashed curves correspond to G-PHD filter and red solid curves correspond to 
G-CPHD filter, (b): True target tracks (solid red) and estimated target tracks 
(dashed with markers) obtained using the G-CPHD filter when Cu, = 2 and 

Pd = 0.9. 


VH. Conclusions 

In this paper we address the problem of multitarget tracking 
using multiple sensors. Many of the existing approaches do 
not make complete use of the multisensor information or are 
computationally infeasible. The contribution of this work is 
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twofold. As our first contribution we derive update equations 
for the general multisensor CPHD filter. These update equa¬ 
tions, similar to the general multisensor PHD hlter update 
equations, are combinatorial in nature and hence computation¬ 
ally intractable. Our second contribution is in developing an 
approximate greedy implementation of the general multisensor 
CPHD and PHD filters based on a Gaussian mixture model. 
The algorithm avoids any combinatorial calculations without 
sacrificing tracking accuracy. The algorithm is also scalable 
since the computational requirements grow linearly in the 
number of sensors as observed from the simulations. 
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Appendix A 

Recursive expression eor partitions 

Let be the collection of all permissible partitions 

of the set {\ < I < s) where partitions are as de- 

hned in the equations ([n])-([T4l). Since the V component 
of a partition is unique given the W components, we do 
not explicitly specify the V component in the recursive 
expression. Let P e be any partition of which 

is given as P = {Wi,W 2 , ■ ■ ■, ,W\p\_i,V}. Let = 

Then we can express using 

p(') and Z'+i as given by the following relation 
■p(ui) _ 

mj+i inin( 7 ni+i,|P|-l) 

u u u u u u 

Pg-ptO ni =0 n 2=0 /2S[l,mi+i] Js[l,|P|-l] 

\Il\=ni 1/21=11,2 l■/|=|/2| 

/in/ 2=0 

U {^3(12)7^12 }i26/2 

BiB{l 2 ,J) 

(39) 

where is the collection of all possible matchings]^ 

from set I 2 to set J. The above relation mathematically 
expresses the fact that for each partition P e and given 
Zj,^\, a new partition belonging to can be constructed 

by adding some new singleton measurement subsets from 
(i.e. extending some existing subsets in 

P by appending them with measurements from (i.e. 

{^B(i2)'“’^i2^}*26/2) and retaining some existing measurement 
subsets (i.e. {Wj}j^j). By this dehnition we have V = 

As a special case for Z = 1, 

mi 

= u (40) 

n =0 

\I\=n 


J) is the collection of all possible one-to-one mappings from set I 

to set J. 


Appendix B 

Functional derivatives 

We now review the notion of functional derivatives which 
play an important role in the derivation of filter update 
equations. A brief background is provided that is necessary 
for the derivations in this paper; for additional details see Q, 
0ch. 11], fig. Let T denote the set of mappings from y to 
K and let A be a functional mapping elements of T to K. Let 
u{jf) and g(y) be functions in T. For the functional A[u], its 
Gateaux derivative along the direction of the function g{y) is 
dehned as © 


dA def,. A[u + e-g] - A[u] 

— [m] = hm-. 

og e 


(41) 


In this paper we are interested in the Gateaux derivatives 
when the function g{y) is the Dirac delta function (5yi(y) 
localized at yi and the corresponding Gateaux derivatives are 
called functional derivatives 0,0 . In this case, the functional 
derivative is commonly written. 


dA SA 

M = 


dS. 


(42) 


-yi '^yi 

If the functional A[u\ is of the form A[u] = f u{y)g{y)dy 
then we have j^[u] = g{yi). 

We can also dehne higher order functional derivatives of 
A[m]. For a set y = {yi,y 2 ,'")yra}, the order derivative 
is denoted by 

(5"A 


5'^ A r -1 def 

5Y ^ SyiSy2...Sy„ 


M- 


(43) 


We call the functional derivative of A[u] with respect 


ZlAi 

SY y 

to the set Y. 

For functionals Ai[u], A 2 [m] and A 3 [u], the product rule 
for functional derivatives ||^ Ch. 11] gives 


— {Ai[w] A2[u] A3[u]} 


= E E 

YiZY Y 2 ZY 

YinY2=0 


5Ai SA 2 . , 5A 3 


Y| |Y| 

= E E E E 

ni=0n2=0 YiSY Y 2 SY 

|Yi|=ni |Y2|=n2 

YinY2=0 


5Ai 6 A 2 ^ 


[“] 


M 3 


■^ 1 -^ 2 ) 


(44) 


M- 


(45) 


As a special case, for the product of two functionals we have 


1^1 


Ml 


{A,[u]A2[u]}=Z E 


M 2 


SY 


n =0 YiSY 

|YiH 


Ml ^ ^SiY-Yi) 


[m] . (46) 


Appendix C 

Probability generating eunctional 

Let u(x) be a function with the mapping u '■ X ^ [0, 1], 
where X is the single target state space. For a set 2f £ A, 
dehne = Let 21 be a random hnite set with 

elements in X and let fs{X) be its probability density 
function. The probability generating functional (PGFL Q) of 
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the random finite set S is defined as the following integral 
transform 

Gs[u] = J f3{X)SX ( 47 ) 

where the integration is a set integral 0- 

For a constant a € [0,1] denote by m ( x ) e a the constant 
function u(x) = a, for all x 6 A”. Let A[a] denote the value 
of the functional A[u\ evaluated at u(x) e a. Recall that for 
a random variable its moments are related to the derivatives 
of its PGR Similarly, the first moment or the PHD function 
of a random finite set is related to the functional derivative 
of its PGFL. For the random finite set S, its PHD function 
i7s(x) is related to the functional derivative of the PGFL 0 
as follows 

SC- 

^ h ( x ) = ^[ 1 ]. ( 48 ) 

The PGF M 3 (t) of the cardinality distribution of the 
random finite set S is obtained by substituting the constant 
function u(x) e t, in the PGFL i.e., M 3 (t) = Gs[t]. For an 
IIDC random finite set with spatial density function C(x) we 
have the relation Gh[m] = Mh((C,u))- 

Appendix D 

Multitarget Bayes eilter 

Let fk+i\k{X\Zl;^) and fk+i\k+i{X\Zl;^^-^) be the predicted 
and posterior multitarget state distributions at time k + 1 and 
let Lk+ij{Z^f_^^\X) denote the multitarget likelihood function 
for the sensor at time k+1. Since the sensor observations 
are independent conditional on the multitarget state, the update 
equation for the multitarget Bayes filter Q is given by 

fk.,\k.i{X\zl;^k.i) - fk.i\kiX\Zl:^k) fl Lk.iAZiJX). 

( 49 ) 

We now define a multivariate functional which is the integral 
transform of the quantity in the right hand side of the above 
equation. Under the conditions of Assumption 1, we can obtain 
a closed form expression for this multivariate functional, which 
on differentiation gives the PGFL of the posterior multitarget 
state distribution. 

Let gj{z),j = l,2,...,s be functions that map the space 
to [0,1] where Z-' is the space of observations of sensor 
j. The intermediate functions gj{z) will be used to define 
functionals and later set to zero to obtain the PGFL of the 
posterior multitarget distribution. Let u(x) be a function 
mapping the state space X to [0,1]. For brevity, denote 
the vector of functions [51,52, as 5i;s and define 

5 j(z) where Z^ £ ZK We define the multivariate 
functional ^[51,52 ,... , 5 s,m] as the following integral trans¬ 
form 

F\g^..,,A = ^ |n/:fc,i.,[5,|X]j jk^^A^\Z\A)bX 

(50) 

where Ck^^A9,\X] = f gf Lfc+i., |X)<5ZL (51) 


Later we will relate the PGFL of the posterior multitarget 
distribution to the derivatives of the functional F[gi-s,u] with 
respect to the sensor observations Z^i- Recall that Cj{z) 
denotes the clutter spatial distribution and Gj{t) denotes the 
PGF of the clutter cardinality distribution for the 5*^ sensor. 
Under Assumption 1 it can be shown that Q 

Ck^iAgj\X] = f gf Lk^iAZ^\X)SZ^ (52) 

= G,((c„5,))4, (53) 

where (x) = 1 - p^x) + AdM Pa, (^) (54) 

Pa, (^) = / Po (z)/ii(z|x)dz. (55) 

Let Gk+i\k[u\ denote the PGFL of the predicted multitarget 
distribution. Using the above relations in ( |50l l we have 

F[g,..sM = f I fk.i\k{X\ZAk)5X 


Since both u and (j)g. are functions defined over the space X, 
we can combine the product of and nj=i write 

{u nj=i 4‘gj)^■ Hence we have 


F[ 9 i-.s,u\ 

= (fl fl j fk.i\k{X\Z^A) 5X 

(57) 


' s \ s 

Y\Gg{{cg,gg))\G[uY\K^ ( 58 ) 

5=1 / J =1 

' S \ S 

Y\GA{03.9a))\M{{T,uY\K))- ( 59 ) 

5=1 / 1=1 


The last two steps result from the definition of the PGFL 
and the assumption that the predicted multitarget distribution 
fk^i\k{X\ZAk) is IIDC. 

Let Gk+i\k+i\u] be the PGFL of the multitarget density 
/fe+i|fc+i(-’^|^ufc+i), and let Dk+i\k+i{x) be the posterior PHD 
function. From 0,1113 we have the following relation 

5F , 

^^TTr[ 0 , 0 ,..., 0 ,M] 

Gk+i\k+i[u] = AA -• (60) 


Since the PHD is the functional derivative of the PGFL, 
from ( |48l l 


Dk+l\k+l{x-) 


6F 


SxSZl'^^ 


[ 0 , 0 ,..., 0 , 1 ] 


SF 


[ 0 , 0 ,..., 0 , 1 ] 


(61) 


Note that the differentiation —J— is with respect to the 

function variable gj and the differentiation ^ is with respect 
to the function variable u. The general multisensor CPHD 
filter update equation is derived by evaluating the functional 
derivatives of F[5i;s,u] in ([60ll and (|6T]). 
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We now define a quantity T and the functionals 
and (f-wYgi-.a^u]- The functional derivatives of can 

be expressed in terms of these quantities. Let 

\ 

(62) 


r=n n c,(z) 




Vj=i 

For W e W, let 

‘Pw[9i-.s,u] 


j=i 


y'r(x)w(x)| n j I n 


\(i,l)iTw 




n Ci(zi) 

{i,l)^Tw 


6F 




where. 


r<'’sn n «,{-). 




I 




and for W € >V, 

g^w[gi:s,u] 

y'r(x)u(x)| n I I n <?!>3,(x)|rfx 


^(z,/)eTvtp 




n c.(zi) 


(69) 


(63) 


(64) 


With these definitions we can prove, via mathematical induc¬ 
tion, the following lemma. 

Lemma 1. Under the conditions of Assumption 1, the func¬ 
tional derivative of F'igi-s,u\ with respect to the multisensor 
observation set is given by 

SF 


Mathematical induction: case I = 1 

We first establish the induction result for the base case, i.e. 
I = 1. Ignoring the time index let the observation set gathered 
by sensor 1 at time fc + 1 be = {z}, z^,..., We 

have, for the case of s sensors 

F[gi-.a,u] = (nQ((cj-,gi)) I M{{r,ufl(j)g^)). (70) 

V=i / 3=1 

Differentiating the above expression with respect to the set 
ZLi we get 
6F 

- ^{(nQCiq.S,))) M((r,«n.#„))) (71) 


— [9i-.8,u] = r Y, 4'p[ffi:s,M] n ‘Pw[gi-.s,u] (65) 

PeP WsP 

where F, T'p[pi:s,'u] and V^w[ 5 i:s,'u] are as defined in ( |62| i, 
and 

Lemma is proved in Appendix 

Appendix E 
Proof of Lemma[T] 

Proof The derivation is based on the approach used by 
Mahler Q to derive multisensor PHD filter equations for the 
two sensor case and its extension by Delande et al. pO) for 
the general case of s sensors. 

Mathematical induction 

We prove using mathematical induction on 1 < Z < s the 
following result. 


3 = 1 


(72) 


since the differential — only differentiates the variable gi. 

If / £ we can express Y £ Z^i as F = {z^ : i € I] 

for some /. We also have Z^i - Y = {z^ : f ^ /}. Using the 
product rule for functional derivatives from (|46]l we have 


I C'i((ci,5i)) M((r, u fl )) 


Pspto WiP 

( 66 ) 


fc+i 

mi c ® A 

= E E n ^((mri</>9i)) .. n 67i((ci,gi)). 

n=0 7c[l,„ii] nZi'lie/ J = 1 <5{zj}i^7 

\I\=n 

(73) 

Now we consider the derivatives of each of the individual 
terms in the above expression. By application of the chain 
rule for functional derivatives ||^ Ch. 11] 

16/ 3 = 1 




(67) 


X n Q((c„5,)) (68) 

\j=i+i / j=i 


= n(AMPd^i(zi) (74) 

3 = 1 *6/ j=2 

In the expression above and in later expressions we have used 
the convention that whenever |/| = 0, Hie/O = 1- Applying 
the chain rule to the second derivative 

—4^Ci((ci,5i)) = C^"^-"\(ci,pi)) nci(z!) (75) 

= (76) 

llci(zj 
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As 

0 , 


before we have used the convention His/O = 1 when |/| = 


IS betore we nave used me convention = i wl 

Thus the right hand side of 0 can be expressed 

mi I S 

n=0/c[l,mil I i=l 

\I\=n 


as 




X7b+1 , 

^^k+l [ Pe'P(f’) 


W^P 

b 


■ri‘>| n C,((c„a,})\ E [Yldp™ 


(('"i ’ 9 3 )) 


\j=b+2 


_P6'P(<>) \f = l 




j=i 


n ‘^w[9v.s,u\\. 

WiP 


expression above we have 

5 




3-]C'b+i((cb+i,p6+i))M(*^l 


j=i 


iil CifZ-j j 

In the double summation above, each set I maps to a partition 
P of the form P = IJie/ {zj} in and vice versa. Hence 
using result in equation ( |40l i of Appendix we have 

'^‘^k+1 Pipm [ 

\i=2 / 3 = 1 WiP ) 

(78) 

Note that for the empty partition P = {V}, there are no 
elements of the form W in P. Hence in the above expression 
we use the convention HivepO = 1 whenever P = {V}. 
Further grouping of the terms gives the compact expression 

"^fe+1 PsPW WsP 

(79) 

Hence the result is established for the case I = 1. 

Mathematical induction: case I = b>l 

Now assuming that the result is true for some / = & > 1, we 
establish that the result holds for / = 5 + 1 < s. Let = 
Z 2 '"^ • • ■ 

SF . , S ( 6F . J 

( 80 ) 

Substituting the result for the case I = b we get 

SF 

[gi-.s,u\ 


n 'fw[gi-.s,u]\ 

WiP j 

mi,+i min(mi,+i,|P|-l) 

= E E E E 

ni=0 712=0 /ic[[l,mb+i| / 2 ^[[l,'mb+i| 

|/i|=ni \l 2 \=' 7 i 2 ‘, /in/2=0 






n g^w[gv.s,u\\ 

\WiP 1 


H'^i2^}i2'^l2 \WiP 

5 


5{'^f^]^iPuI2 


Ch+i{{cb+i,gb+i))\- (83) 


The second summation above is restricted to the 
limit min(mb+i,|P| - 1) because the derivatives of 
Y\wiP ‘fw^gi-.sju] for 712 > |7^l“ 1 ^6 zero. Now considering 
each of the individual derivatives above we have 

<5 




"=i 

(84) 


n=l 


Denote P = {VFi, W2 ,..., W|p|_i, L} for notational conve¬ 
nience. Then we have 


n ‘pw[9i-.s,u] 


r(&) ^ [ (81) 


1*26/2 \WiP 


E E 1 X 

Js[l.|P|-l]B6t3(/2.J) ( Vjs'J / 

NHL] 

I n (85) 

where B{l 2 , J) is the collection of all possible matchings from 
set I 2 to set J and we define the measurement subset = 
VFpi,. * u zj+i. Also 


^(■^2) ‘-‘12 

s 


{Cb+i{{cb+i, gb+i))) 


(82) 


'^{^i^^}*^/lU/2 

= n C5.i(zri). (86) 

iS‘/lU/2 

Combining the three derivatives into the right hand side of 
expression ([8^ we get 


Let Ii £ [l,mb+i] and I 2 £ [1,to;,+i] such that /in/2 = 0. 
Then we can express Yi £ and L2 - ^kti satisfying 

Yi n 1X2 = 0 as Yi = {z^+^ : i e Bj and Y 2 = {z^^^ : i e 
12} respectively. Applying the product rule from (|45]l to the 
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n C,{{c,,g,))\{ n 

\j=b+2 




rrib+i min(mi,+ i,|P|-l) 

I I I I I I 

p^'pib) ni =0 n 2=0 /ic[[l,m 5 +i| /2^[[l,'/^&+i| Je|Il,|P|-l| 

|/i|=ni \l2\=n2', /in/2=0 |J|=|/2l 


L(mi +1 m "2)((cb+i,gf,+ l))(nV5lV,[5l:s,M] )x 

,J) I \jtJ } 


I 


n V5{z‘'+i}[5i:s>'“] n ^wi>[9i--s,u\ 

_ ^ r_ / \ ^- ^ r_ ^2 


Uie/i 


22^/2 






P<=p(<>+i) [ \j=i y 

fl p((9>5j)) I </>g„))x 


\i=h+2 


n ‘/5W-[51 :s,m] 1 

VKeP / 


VKseP 


■^fe+1 


W^sP 


i5x 


r Y, ^'pM n ‘•pw[u\ 


= rE me: 

PeP \j=l 


PeP H"eP 

(^Tij-|P|i) 


(0) 


j=l VKeP 


Applying the product rule for set derivatives from < |46| ) 


5x 


Afd^l ^H{r,uflqd)) n 


Jx 


j = l WeP 


^^((r,unq^)) i n mM + 


j=l I WeP 


n (93) 


j = l (WsP 

Evaluating the individual derivatives above and substituting 
the constant function u(x) e 1, we get 


(87) 


Using result of Appendix we can simplify the multiple 
summation term and write 


— |m( 1 I ^((r,wng^))| =M(l^9('^)^(-x)fl(?^(x) 

I J=i J U.1 J=i 

(94) 

^jn =(n <^w]i Y ^(x)pw(x)') (95) 

ox [^r^p \WiP / VtVeP / 

where 7, dw and pvv(x) are dehned in ( [T7] i, ( [T9] ) and ( |20| ) 
respectively. We note that for the empty partition P = {V} 
there are no elements of the form W in P. In this case the 
derivative in equation ( |95l l is zero since the quantity being 
differentiated is a constant equal to 1. To have a compact 
representation of the update equations we use the convention 
T,WipO - 0 when P = {V}. Hence we have 


= E n ^w[9i:s,ul (88) 

Pg-pCb+l) 


SF 


SxdZ 


^[o,o,...,o,i] = r E 


(mj-|P|j) 


(0) 


fc+1 


PeP\g=l 


Hence we have established the result stated in ( |66| l using the 
method of mathematical induction. We obtain the result of 
Lemma [T] by substituting I = s in this result. 

□ 


Appendix F 
Proof of Theorem[T] 

For brevity denote T'p[0,0,..., 0, u] = 'I'p[u] and 

(/?w[0,0,..., 0 ,m] = ipw[u]- Substituting gj = 0 for j = 
1,2,... ,s in the result of Lemma we get 

SF 

-^[0,0,...,0,u] = r E 'I'pM Yl^w[ul (89) 

^ Pep Wc p 


|m(1^9(^)^(x) ng^(x) [ n ^(1^] + 

[ j=l \WiP / 

^(l^l-^)(7)( n dw\i Y ^(x)pvf(x)) I (96) 

VtVeP / VVFeP / j 

= r E («:pM(l^9 Yl dw] [r-(x) n(?^(x)'| 

PeP \ WiP / \ j=l j 

+ r E (ppMd^l-i) n ( E Kx)ptv(x)) (97) 

PeP \ WiP / VtVeP / 

where Kp is dehned in ( [T^ . Substituting u(x) e 1 in equation 
we have 


SF 


^[ 0 , 0 ,..., 0 ,l]=r E n dw- (98) 


PHD update 

Differentiating equation (89i with respect to set {x} we 
have 

SxSZY\ [0> 0, ■ • •, 0, u] = — I[0,0,..., 0, u]| (90) 


d^k+l PeP WeP 

Dividing (j^ by ( |98| ) and using the dehnition of PHD 
from ( jhTj l, we get 

SF , 

^[ 0 , 0 ,..., 0 , 1 ] 

- (99) 


, ^ SxSZP 

i9fc+i|fc+i(x) = — 


(91) 




— |m(I^I i(((r,ung^)) n V’wlu]}- (92) 


d^k+l 

= r(x) I ao fl 9d(x) + E ( E Pif(x)] [ (100) 

[ j = l PeP VtVeP / j 

where ao and ap are as given in ( |2T] i and ( |22] l. 

Cardinality update 

We now derive the update equation for the posterior car- 
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dinality distribution. Using the expression for the posterior 
probability generating functional in ( |60l ) and the results of ( |8^ 
and ([98]l we have 


Gk+l\k+l [l*] 


z n '^w[u\ 

PeP WeP 


^ Yl d^- 

PeP W^P 


( 101 ) 


We thus have 


Pk+l\k+l (^) 

Pfc+l|fc (^) 


E («p 

PeP \ 


PeP 

|P|<n+l 


(n-|P| + l)! 


n-|P|+l 


n 

WiP / 


^ n dw 


PiV 


WiP 


( 110 ) 


The probability generating function of the pos¬ 

terior cardinality distribution is obtained by substituting the 
constant function m(x) e f in the expression for Gfc+i|fc+i[w]. 
Thus 


d^k+l\k+l{t) - 


z n pw[t] 

PeP WeP 


^ ppM0^1-i) n dw' 

PeP WeP 


For constant t we have 

n '^w[t ]=n dw 


( 102 ) 


(103) 


WiP 


WiP 


and T'p[f] = ( fj i (104) 


Vi=i 


Since is the PGF corresponding to the cardinality 

distribution Pk+i\k+i{n), 




1 

n! 


E n <1 

a PeP 


w 


WeP 


(105) 

(106) 


i=0 


' dt^ Y. n dw 

PiP WiP 

Pep \j=l / WeP J 


n\ Y n dw 


PeP 


WeP 


(107) 


Evaluating the derivative we get 


dt^ ^ 


t=o 


.{n-\P\ + l)l 


if p < |P| - 1 
M(”)(0)7"-l^i+i ifp>|P|-l. 


(108) 


We also have M1”1(0) = n!pfc+i|fc(n), hence 


z 


ptp (n-|P| + l)! ^ 

|P|<n+l 


^n-|PFl Yl dw 


W^P 


Y kpMO^I-1) Yl dw 


PeP 


WiP 


( 109 ) 


where up is as defined in ( [TSl ). 
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